Project Question 2

0 - Setup

if (!require("pacman")) 
    install.packages("pacman")

# Use pacman::p_load to install and load CRAN packages
pacman::p_load(
    countdown,
    # datasauRus,
    dplyr,
    elevatr,
    gganimate,
    ggplot2,
    ggthemes,
    ggspatial,
    gifski,
    glue,
    grid,
    gridExtra,
    here,
    hms,
    knitr,
    lubridate,
    metR,
    openintro,
    readr,
    rnaturalearth,
    rnaturalearthdata,
    rnaturalearthhires,
    terra,
    scales,
    sf, 
    tidyverse,
    viridis
)

# Handle GitHub package separately
if (!require("dsbox")) {
      # Install devtools if not present
      if (!require("devtools")) 
          install.packages("devtools")
devtools::install_github("tidyverse/dsbox")
      library(dsbox)
}

ggplot2::theme_set(ggplot2::theme_light(base_size = 12))
options(width = 65)

knitr::opts_chunk$set(
  fig.width = 7, # 7" width
  fig.asp = 0.618, # the golden ratio
  fig.retina = 3, # dpi multiplier for displaying HTML output on retina
  fig.align = "center", # center align figures
  dpi = 300 # higher dpi, sharper image
)

save_figs <- TRUE

vesuvius_data <- read.csv(here("data","vesuvius.csv"))

vesuvius_data_filtered <- vesuvius_data |>
    drop_na(
        latitude, 
        longitude, 
        depth_km, 
        duration_magnitude_md
    ) |>
    mutate(
        depth_abs = abs(depth_km),
        depth_boundary_1km = case_when(
            depth_km < 1 ~ "< 1 km",
            depth_km < 2 ~ "< 2 km",
            depth_km < 3 ~ "< 3 km",
            depth_km < 4 ~ "< 4 km",
            depth_km < 5 ~ "< 5 km",
            TRUE ~ "> 5 km"),
        
        depth_boundary_1km = factor(
            depth_boundary_1km, 
            levels = c("< 1 km", 
                       "< 2 km", 
                       "< 3 km", 
                       "< 4 km", 
                       "< 5 km", 
                       "> 5 km")),
        
        depth_boundary_0.5km = case_when(
            depth_km < 0.5 ~ "< 0.5 km",
            depth_km < 1 ~ "0.5 ~ 1 km",
            depth_km < 1.5 ~ "1 ~ 1.5 km",
            depth_km < 2 ~ "1.5 ~ 2 km",
            depth_km < 2.5 ~ "2 ~ 2.5 km",
            depth_km < 3 ~ "2.5 ~ 3 km",
            depth_km < 3.5 ~ "3 ~ 3.5 km",
            depth_km < 4 ~ "3.5 ~ 4 km",
            depth_km < 4.5 ~ "4 ~ 4.5 km",
            depth_km < 5 ~ "4.5 ~ 5 km",
            TRUE ~ "> 5 km"),
        
        depth_boundary_0.5km = factor(
            depth_boundary_0.5km, 
            levels = rev(c("< 0.5 km", 
                       "0.5 ~ 1 km", 
                       "1 ~ 1.5 km", 
                       "1.5 ~ 2 km", 
                       "2 ~ 2.5 km", 
                       "2.5 ~ 3 km",
                       "3 ~ 3.5 km", 
                       "3.5 ~ 4 km", 
                       "4 ~ 4.5 km", 
                       "4.5 ~ 5 km", 
                       "> 5 km")))
    ) |>

    filter(
        !is.na(depth_boundary_0.5km)
    ) |>

    mutate(
        depth_rescaled = rescale(depth_km, to = c(1, 11)),
        mag_rescaled = rescale(duration_magnitude_md, to = c(0.4, 0.9))
    ) |>

    arrange(
        depth_boundary_0.5km, 
        desc(depth_rescaled)
    ) |> glimpse()

Rows: 8,475
Columns: 16
$ event_id              <int> 34136, 5595, 13138, 1499, 21103, …
$ time                  <chr> "2024-04-14T16:05:54Z", "2021-08-…
$ latitude              <dbl> 40.86050, 40.82580, 40.82570, 40.…
$ longitude             <dbl> 14.34917, 14.43120, 14.43000, 14.…
$ depth_km              <dbl> 9.35, 5.11, 4.78, 4.42, 4.00, 4.0…
$ duration_magnitude_md <dbl> 2.30, 1.60, 2.18, 2.50, 1.70, 0.0…
$ md_error              <dbl> 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3…
$ area                  <chr> "Mount Vesuvius", "Mount Vesuvius…
$ type                  <chr> "earthquake", "earthquake", "eart…
$ review_level          <chr> "preliminary", "preliminary", "pr…
$ year                  <int> 2024, 2021, 2021, 2019, 2014, 201…
$ depth_abs             <dbl> 9.35, 5.11, 4.78, 4.42, 4.00, 4.0…
$ depth_boundary_1km    <fct> > 5 km, > 5 km, < 5 km, < 5 km, <…
$ depth_boundary_0.5km  <fct> > 5 km, > 5 km, 4.5 ~ 5 km, 4 ~ 4…
$ depth_rescaled        <dbl> 11.000000, 6.460385, 6.107066, 5.…
$ mag_rescaled          <dbl> 0.8069767, 0.7255814, 0.7930233, …

max_depth <- max(vesuvius_data$depth_km, na.rm = TRUE)
min_depth <- min(vesuvius_data$depth_km, na.rm = TRUE)

print(paste("Maximum depth (km):", max_depth))

[1] "Maximum depth (km): 9.35"

print(paste("Minimum depth (km):", min_depth))

[1] "Minimum depth (km): 0.01"

max_mag <- max(vesuvius_data$duration_magnitude_md, na.rm = TRUE)
min_mag <- min(vesuvius_data$duration_magnitude_md, na.rm = TRUE)

print(paste("Maximum Magnitude:", max_mag))

[1] "Maximum Magnitude: 3.1"

print(paste("Minimum Magnitude:", min_mag))

[1] "Minimum Magnitude: -2"

original <-
    ggplot(
        vesuvius_data_filtered, 
        aes(x = longitude, 
            y = latitude)
    ) +
    
    geom_point(
        aes(color = duration_magnitude_md,
            alpha = depth_abs),
        size = 3
    ) +                     
    
    scale_color_viridis_c(
        option = "inferno",
        name = "Magnitude"
    ) +
    
    scale_alpha_continuous(
        name = "Depth (km)",
        range = c(0.75, 0.25)
    ) +

        scale_x_continuous(
            limits = c(14.34, 14.48), 
            expand = c(0, 0)
        ) +
        scale_y_continuous(
            limits = c(40.78, 40.88), 
            expand = c(0, 0)
        ) +
    
    labs(
        title = "Seismic Events at Mount Vesuvius",
        subtitle = "Point Color by Magnitude\nPoint Transparency by Depth",
        x = "Longitude",
        y = "Latitude",
        caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
    ) +

    theme(
        plot.title.position = "plot",
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.direction = "horizontal",
        legend.box = "vertical",  # stack legends on top of each other
        legend.box.just = "left",
        legend.spacing.x = unit(0.3, "cm"),  # tighter spacing
        legend.spacing.y = unit(0.3, "cm"),  # tighter spacing
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 10),
        legend.background = element_rect(fill = "transparent", colour = NA),
        legend.box.background = element_rect(fill = "transparent", colour = NA)
    ) +

    guides(
        color = guide_colorbar(  # Magnitude
            direction = "horizontal",
            title.position = "top",
            barwidth = unit(2.8, "cm"),
            barheight = unit(0.4, "cm"),
            order = 1),
    alpha = guide_legend(
        direction = "horizontal",
        title.position = "top",
        label.position = "bottom",
        override.aes = list(alpha = c(0.75, 0.5, 0.25)),
        keywidth = unit(0.4, "cm"),
        keyheight = unit(0.4, "cm"),
        order = 2)
    )
        
original

if (save_figs) {
    ggsave("images/original_vesuvius_plot.png", original)
}

histogram <- 
    ggplot(
        vesuvius_data_filtered, 
        aes(x = duration_magnitude_md)
    ) +

    geom_histogram(
        binwidth = 0.2, 
        fill = 'red', 
        color = "white"
    ) +

    facet_wrap(
        ~ depth_boundary_0.5km, 
        scales = "free_y"
    ) +

    labs(
        title = "Histogram of Duration Magnitude",
        subtitle = "Distribution Facets by Depth",
        x = "Magnitude (md)",
        y = "Count",
        caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
    ) +
    
    theme(
        plot.title.position = "plot"
    )

histogram

if (save_figs) {
    ggsave("images/histogram_vesuvius_plot.png", histogram)
}

# Define Vesuvius bounding box (in WGS84 lon/lat)
vesuvius_bbox <- matrix(c(14.33, 40.78, 14.48, 40.89), ncol = 2, byrow = TRUE)

# Convert bbox to sf polygon
vesuvius_sf <- st_as_sf(
  st_sfc(st_polygon(list(rbind(
    c(vesuvius_bbox[1,1], vesuvius_bbox[1,2]),
    c(vesuvius_bbox[2,1], vesuvius_bbox[1,2]),
    c(vesuvius_bbox[2,1], vesuvius_bbox[2,2]),
    c(vesuvius_bbox[1,1], vesuvius_bbox[2,2]),
    c(vesuvius_bbox[1,1], vesuvius_bbox[1,2])
  ))), crs = 4326))

# Get elevation data from AWS Terrain Tiles (zoom 10)
elev_raster <- get_elev_raster(
    locations = vesuvius_sf, 
    z = 10, 
    clip = "locations")

# Convert to dataframe for ggplot
elev_df <- as.data.frame(
    elev_raster, 
    xy = TRUE)

colnames(elev_df) <- c("x", "y", "elevation")

static_plot <- ggplot() +
    metR::geom_contour_fill(
        data = elev_df,
        aes(x = x, 
            y = y,
            z = elevation/1000, 
            group = NULL),
        bins = 10
    ) +

    scale_fill_gradientn(
        colors = gray.colors(10, start = 1.0, end = 0.3),
        name = "Elevation (km)",
        limits = c(-0.111, 1.246),  # Match elevation range from elev_df (km)
        breaks = seq(0, 1.2, by = 0.3)
    ) +

    geom_point(
        data = vesuvius_data_filtered,
        aes(x = longitude, 
            y = latitude,
            color = duration_magnitude_md,
            size = depth_rescaled),
        alpha = 0.75
    ) +

    scale_size_continuous(
        name = "Depth (km)",
        breaks = c(2.5, 5, 7.5),
        range = c(1, 11)  # adjust to control point size
    ) +

    scale_color_viridis_c(
        option = "magma", 
        name = "Magnitude (md)"
    ) +

    scale_x_continuous(
        limits = c(14.34, 14.48), 
        expand = c(0, 0)
    ) +
    scale_y_continuous(
        limits = c(40.78, 40.88), 
        expand = c(0, 0)
    ) +

    labs(
        title = "Seismic Events at Mount Vesuvius",
        subtitle = "Point Color by Magnitude\nPoint Size by Depth",
        x = "Longitude",
        y = "Latitude",
        caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
    ) +

    theme(
        plot.title.position = "plot",
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.direction = "horizontal",
        legend.box = "vertical",  # stack legends on top of each other
        legend.box.just = "left",
        legend.spacing.x = unit(0.3, "cm"),  # tighter spacing
        legend.spacing.y = unit(0.3, "cm"),  # tighter spacing
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 10),
        legend.background = element_rect(fill = "transparent", colour = NA),
        legend.box.background = element_rect(fill = "transparent", colour = NA)
    ) +

    guides(
        fill = guide_colorbar(  # Elevation
            direction = "horizontal",
            title.position = "top",
            barwidth = unit(2.8, "cm"),
            barheight = unit(0.4, "cm"),
            order = 1),
        color = guide_colorbar(  # Magnitude
            direction = "horizontal",
            title.position = "top",
            barwidth = unit(2.8, "cm"),
            barheight = unit(0.4, "cm"),
            order = 2),
        size = guide_legend(  # Depth
            direction = "horizontal",
            title.position = "top",
            label.position = "bottom",
            override.aes = list(alpha = 0.7),
            keywidth = unit(0.15, "cm"),
            keyheight = unit(0.4, "cm"),
            order = 3),
        alpha = "none"
    )

static_plot

if (save_figs) {
    ggsave("images/static_vesuvius_plot.png", static_plot)
}

ggplot2::theme_set(ggplot2::theme_light(base_size = 14))
options(width = 65)

knitr::opts_chunk$set(
  fig.width = 7, # 7" width
  fig.asp = 0.618, # the golden ratio
  fig.retina = 2, # dpi multiplier for displaying HTML output on retina
  fig.align = "center", # center align figures
  dpi = 150 # higher dpi, sharper image
)

vesuvius_anim_plot <- ggplot() +

    geom_point(
        data = vesuvius_data_filtered,
        aes(x = longitude,
            y = latitude,
            color = duration_magnitude_md,
            size = depth_rescaled),
        alpha = 0.75
    ) +

    scale_x_continuous(
        limits = c(14.34, 14.48), 
        expand = c(0, 0)
    ) +
    scale_y_continuous(
        limits = c(40.78, 40.88), 
        expand = c(0, 0)
    ) +

    scale_size_continuous(
        name = "Depth (km)",
        breaks = c(2.5, 5, 7.5),
        range = c(1, 11)
    ) +

    scale_color_viridis_c(
        option = "magma", 
        name = "Magnitude (md)"
    ) +

    labs(
        title = "Seismic Events at Mount Vesuvius",
        subtitle = "Point Color by Magnitude\nPoint Size by Depth: {closest_state}",
        x = "Longitude",
        y = "Latitude",
        caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
    ) +

    theme(
        plot.title.position = "plot",
        plot.title = element_text(margin = margin(b = 23), 
                                  size = 16),
        axis.text = element_text(size = 12),
        legend.direction = "horizontal",
        legend.box = "vertical",  # stack legends on top of each other
        legend.box.just = "left",
        legend.spacing.x = unit(0.3, "cm"),  # tighter spacing
        legend.spacing.y = unit(0.3, "cm"),  # tighter spacing
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 10),
        plot.caption = element_text(margin = margin(t = 30), size = 10, hjust = 1),
        legend.background = element_rect(fill = "transparent", colour = NA),
        legend.box.background = element_rect(fill = "transparent", colour = NA)
    ) +

    guides(
        color = guide_colorbar(
            direction = "horizontal",
            title.position = "top",
            barwidth = unit(2.8, "cm"),
            barheight = unit(0.4, "cm"),
            order = 1),
        size = guide_legend(
            direction = "horizontal",
            title.position = "top",
            label.position = "bottom",
            override.aes = list(alpha = 0.7),
            keywidth = unit(0.15, "cm"),
            keyheight = unit(0.4, "cm"),
            order = 2),
        alpha = "none"
    ) +

    transition_states(
        depth_boundary_0.5km,
        transition_length = 3,
        state_length = 2
    ) +

    shadow_mark(
        past = TRUE, 
        future = FALSE, 
        alpha = 0.2
    ) +

    enter_fade() +
    exit_fade()

gganimate::animate(vesuvius_anim_plot)

if (save_figs) {
    gganimate::anim_save("images/vesuvius_animation.gif", vesuvius_anim_plot)
}